Molecular dynamics simulations of ballistic annihilation 
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Using event-driven molecular dynamics we study one- and two-dimensional ballistic annihilation. 
We estimate exponents ^ and 7 that describe the long-time decay of the number of particles (n(t) 
t~^) and of their typical velocity {v{t) ~ t~^). To a good accuracy our results confirm the scaling 
relation ^-1-7=1. In the two-dimensional case our results are in a good agreement with those 
obtained from the Boltzmann kinetic theory. 
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Kinetics of reacting particle systems is recently inten- 
sively studied. Such systems model a wealth of different 
phenomena and in addition offer an important testing 
ground for nonequilibrium statistical mechanics. In lat- 
tice formulation such systems are much easier to analyze 
and in some cases even exact results could be obtained. 
More realistic, off-lattice (ballistic) versions are at the 
same time much more difficult to examine, but even here 
some analytical insight is possible 

Despite being one of the simplest reacting particle sys- 
tems, single-species ballistic annihilation, A -\- A ^ 0, 
so far evaded exact solution. Consequently, our under- 
standing of the kinetics of this process comes mainly from 
approximate methods. To describe ballistic annihilation 
one usually calculates density of particles n{t) and their 
typical velocity v(t) defined using the second moment of 
the time-dependent velocity distribution f{v,t) 
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v'^f{v, t)dv 



(1) 



Of particular interest are then the exponents f and 7 that 
describe the asymptotic decay of these quantities 



(2) 



A simple power counting or scaling arguments can be 
used to show that 



e + 7 = l 



(3) 



Although a rigorous justification of this relation is miss- 
ing, Piasecki et al. [3 analysing BBGKY-like hierarchy 
presented convincing arguments supporting relation 
A lot of efforts was made to calculate f and 7 |2i f 
lEIEQ' However, it would be desirable to obtain more 
accurate estimations of these exponents since it could 
increase our understanding in this field. For example, 
some arguments suggest '5| that in dimensions d > 1 
Boltzmann kinetic theory provides an accurate descrip- 
tion of such a process. However, a direct verification of 
this hypothesis, using molecular dynamics simulations, is 
not yet fully completed mainly due to too small number 
of particles that were taken into account. Let us also no- 
tice that current estimations of f for d = 1,2 and 3 are 



surprisingly close to the result obtained in the so-called 
reaction-controlled hmit = 4d/(4d + 1)) Q. Further 
check whether this exponent could be given as such a 
simple fraction would be desirable. Let us notice that 
a closely related process, namely ballistic aggregation, is 
solvable in d = 1 and some of its exponents are indeed 
simple fractions 

In the present paper we describe results of extensive 
molecular dynamics simulations of ballistic annihilation 
in d — 1 and 2. Since in such a process particles interact 
(i.e., annihilate) only upon collision efficient simulations 
can be made with the use of the so-called event-driven 
dynamics One of the difficulties in this technique is 
the search for the time of the nearest collision. Most ef- 
ficient algorithms to locate such an event arrange data 
on the heap tree Il0|. For event-driven simulations for 
d > 1 systems it is also essential to include sectorization 
and search for a collision partner only within a given sec- 
tor and its nearest neighbours. Heap tree searching and 
sectorization were already used in event-driven dynamics 
of various hard-core or granular systems |TlL Il2l hj but 
to our knowledge not in the ballistic annihilation. Since 
it substantially improves the numerical performance, our 
method implements these techniques. The fact that at 
the collision particles annihilate, simplifies the algorithm 
comparing to the non-annihilating systems. This is be- 
cause there is no need to examine post-collision events of 
colliding particles. In the following we discuss obtained 
numerical results. 

In the d — 1 case we place N point-like particles on 
a line interval of length L = N. Such a choice keeps 
the initial density of particles n — N/L constant (and 
equal to unity) and that will allow us the comparison of 
results for different values of N. Initially velocities have 
a distribution f{v,t = 0) that either will be gaussian 
or it will have a v = singularity f{v,t = 0) ~ 
with a characteristic exponent fj.. Periodic boundary con- 
ditions are used. We measured the density of particles 
n{t) at time t and their typical velocity v(t) defined using 
Eq. (Q. To estimate the exponents 1^ and 7 the system 
must reach the asymptotic regime. However, for finite 
N the long-time behaviour is affected by poor statistics 
as well as by finite-size effects and that hinders precise 
determination of ^ and 7. To take these effects into ac- 
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FIG. 1: The exponent ^ as a function of t for TV = 10^, 5 ■ 
10®,10'^,2 ■ 10^, and 3 • lO'' (thick line). Initially, velocities 
have a gaussian distribution (/x = 0). 

count we calculated time-dependent exponents ^(t) and 
7(t) using results (i.e., n(t) and v{t)) spreading over one 
decade around a given value of t. 

For the gaussian initial velocity distribution our re- 
sults are shown in Fig. I1I2I Presented values are aver- 
ages over 3 • 10* (for N = 3- 10^) to lO'^ (for N = 10<^) 
samples. Although examined systems were rather large 
{Nmax — 3 • 10^), it is still difficult to make precise 
evaluation of exponents. Our results, ^ — 0.805(2) and 
7 = 0.195(2), satisfy Eq. ©, but we admit that the esti- 
mation of errors is based mainly on the visual inspection 
of data and the hope that further increase of the system 
size will not change much the final estimations. Although 
the difference is very small, such a value of ^ seems to 
exclude the possibility ^ — 0.8 obtained (exactly) for 
the d = 1 annihilation process in the reaction-controlled 
limit. Our estimation of ^ also slightly differs from previ- 
ous molecular dyiiamics simulations result 0.785(5) made 
by Rey et al. |l5| albeit with much smaller number of 
particles (they simulated systems with up to 262144 par- 
ticles). Let us notice, that computing time in the search 
algorithm of Rey et al. increases as 0{N^/Hn{N)) while 
for the heap tree search (that we used) it increases only as 
0{Nln{Nj). Various computations based on the solution 
of the Boltzmann equation, and thus based on the molec- 
ular chaos hypothesis, predicts ^ close to 0.77 ^i^]. Dis- 
agreement between this result and our estimation shows 
that in d = 1 Boltzmann kinetic theory has only limited 
validity. 

To check for a possible multi-scaling we calculated the 
typical velocity v'{t) using the fourth-moment analogue 
of Eq. 1^. Fig.|3]shows the exponent 7' that describes the 
decay of v'{t). Within our numerical precision 7 = 7' and 
thus v{t) as defined in Eq. is the only characteristic 
velocity in this problem. 

We also performed calculations for the singular initial 
velocity distribution with fi = 1/2. In this case we esti- 
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FIG. 2: The exponent 7 as a function of t iox N = 10®, 5 • 
10^,10^,2 • 10'^, and 3 • lO'^ (thick line). Initially, velocities 
have a gaussian distribution (/i = 0). 
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FIG. 3: The time dependence of the exponent 7' that de- 
scribes the decay of the typical velocity defined via the fourth 
moment. Calculations were made for = 10^ 5 • 10*', 10^ 2 ■ 
10^, and 3T0^ (thick line). Initially, velocities have a gaussian 
distribution (/i = 0). 



mate (see Figs.|3El) i = 0.585(2) and 7 = 0.415(2). This 
result satisfies Eq. (PI and improves over previous esti- 
mations of_£ reported in the literature that range from 
0.5 to 0.6 lIUlIl. 

In the one dimensional case topological constraints im- 
ply that for a given particle potential partners for a col- 
lision are only two of its nearest neighbours. In two di- 
mensions this is no longer the case and in principle any 
pair of particles can collide. For large TV it leads to a pro- 
hibitively large number of potential collisions that should 
be examined but only very few of them will actually take 
place. To overcome this difficulty one can divide available 
space into sectors and look for potential collisions only 
within a sector or neighbouring sectors. More details on 
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FIG. 4: The exponent ^ as a function of t for iV = 10*^, 3-10'', 
and 10^ (thick line). Initially, velocities have a distribution 
with fj. — 1/2 singularity. 
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FIG. 6: The exponent ^ as a function of t for A'^ = 
300^ 500^ 1000^ 2000^ and 3000^ (thick line). Initially, 
velocities have a gaussian distribution (jj, = 0). 
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FIG. 5: The exponent 7 as a function of f for iV = 10^ , 3 • 10*^ , 
and 10'^ (thick line). Initially, velocities have a distribution 
with fi = 1/2 singularity. 
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FIG. 7: The exponent 7 as a function of t for = 
300^ 500^, 1000^ and 3000^ (thick line). Initially, veloci- 
ties have a gaussian distribution (/i = 0). 



this technique can be found elsewhere 0, . 

For the gaussian initial velocity distribution results of 
our simulations are shown in Figs l6l7l In two dimensions 
size of particles r relative to the linear system size L be- 
comes a relevant variable that might affect e.g., the time 
needed to reach asymptotic regime. Presented results are 
obtained for the packing fraction / — irr^N/Lp' = 0.0079 
but similar results were obtained for / = 0.031. After av- 
eraging over 10"^ — 10^ samples, we estimate ^ = 0.872(2) 
and 7 = 0.129(2). Such results satisfy Eq. Qand are 
in a very good agreement with calculations based on the 
Boltzmann equation 0, Q, confirming thus validity of 
molecular chaos hypothesis in this case. Let us notice 
that our estimation of ^ definitely excludes the reaction- 
controlled value ^ = 8/9. Moreover, the maximal number 
of particles examined in our approach (9 • 10^) is almost 



20 times larger that for previously reported d = 2 event- 
driven simulations j^. 

We also made calculations for the /i = 1/2 case and 
the results are presented in Fig. I8I9I We estimate ^ = 
0.83(1) and 7 — 0.17(1). Similar value of ^ was recently 
obtained using numerical integration of the Boltzmann 
equation % In Table CI we collect all our final results. 

In conclusion, implementing a heap-tree search algo- 
rithm and sectorization we made extensive event-driven 
molecular dynamics simulations of the ballistic annihila- 
tion. Obtained estimations of decay exponents ^ and 7 
obey the scaling relation ^ + 7 = 1 and improve over pre- 
viously reported results. Our calculations suggest that 
the typical velocity as defined via the second moment 
sets the only characteristic scale in the problem. Very 
good agreement of our simulations with predictions of 
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TABLE I; Exponents ^ and 7 as calculated in the present 
paper. Estimations of ^, based on the analysis of the Boltz- 
mann equation [3.ll4| and in the reaction-controlled limit Q 
are presented in the last three columns. 



FIG. 8: The exponent ^ as a function of t for A'^ — 
300^ 500^ 1000^ 2000^ and 3000^ (thick line). Initially, 
velocities have a distribution with /i = 1/2 singularity. 
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Boltzmann kinetic theories in the d = 2 case confirms 
validity of the latter approach. In the d = 1 case small 
differences between results of these two approaches ex- 
ist. It would be interesting to extend our approach to 
some other model phenomena whose kinetics is still not 
fully understood as e.g., ballistic aggregation [l. or 
probabilistic ballistic annihilation [l7j . 
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FIG. 9: The exponent 7 as a function of t for N — 
500^ 1000^ 2000^, and 3000^ (thick line). Initially, velocities 
have a distribution with fi = 1/2 singularity. 
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